function [J B]=Q4_cal_B(r,s,xy)
Nrs=0.25*[-(1-s) 1-s 1+s -(1+s);
         -(1-r) -(1+r) 1+r 1-r];
J = Nrs * xy;
A = 1/det(J)*[J(2,2) -J(1,2) 0 0
              0 0 -J(2,1) J(1,1)
              -J(2,1) J(1,1) J(2,2) -J(1,2)];
N1r=-0.25*(1-s);N2r=0.25*(1-s);N3r=0.25*(1+s);N4r=-0.25*(1+s);
N1s=-0.25*(1-r);N2s=-0.25*(1+r);N3s=0.25*(1+r);N4s=0.25*(1-r);
G = [N1r 0 N2r 0 N3r 0 N4r 0
     N1s 0 N2s 0 N3s 0 N4s 0
     0 N1r 0 N2r 0 N3r 0 N4r
     0 N1s 0 N2s 0 N3s 0 N4s];
B = A * G;
end